proc format;
value trtfmt 1='Control' 2='Added Water';
/* I need Order for plotting later on--it could
be created 'on the run' as well */
data a;
input Trt Run Order Voltage Current;
TrtCell=10*Trt+Run;
label trtcell="Treatment/Cell";
format Trt trtfmt.;
datalines;
1 1 1 .60 20.3
1 1 2 .45 32.3
1 1 3 .50 29.2
1 1 4 .65 15.9
1 1 5 .55 25.3
1 1 6 .70 10.9
1 1 7 .75  6.59
1 1 8 .80  2.53
1 2 1 .60 21.2
1 2 2 .45 32.8
1 2 3 .50 29.8
1 2 4 .65 16.4
1 2 5 .55 26.0
1 2 6 .70 11.1
1 2 7 .75  6.66
1 2 8 .80  2.53
1 3 1 .60 21.4
1 3 2 .45 33.3
1 3 3 .50 30.1
1 3 4 .65 16.5
1 3 5 .55 26.2
1 3 6 .70 11.2
1 3 7 .75  6.66
1 3 8 .80  2.43
2 1 1 .60 19.5
2 1 2 .45 31.6
2 1 3 .50 28.6
2 1 4 .65 15.6
2 1 5 .55 25.0
2 1 6 .70 10.6
2 1 7 .75  6.36
2 1 8 .80  2.23
2 2 1 .60 19.7
2 2 2 .45 31.4
2 2 3 .50 28.5
2 2 4 .65 15.7
2 2 5 .55 24.9
2 2 6 .70 10.6
2 2 7 .75  6.35
2 2 8  .80  2.20
;
run;
*Plot the data;
proc sort data=a; by trt run voltage;

proc sgplot data=a;
series x=voltage y=current/group=trtcell;
run;

proc sgpanel data=a;
panelby trt / rows=1 columns=2;
series x=voltage y=current/group=run;
run;

/* Use GLM to produce residuals */
ods graphics on;
proc glm data=a plots=all;
class Trt Run Voltage;
model Current=Trt Run(Trt) Voltage Trt*Voltage;
random Run(Trt)/test;
output out=outa r=resid; run;
ods graphics;

/* Plot the ordered residuals */
proc sort data=outa; by Trt Run order; proc print data=outa; run;
proc gplot data=outa; by Trt Run;
plot resid*order=1/haxis=axis2 vaxis=axis1 vref=0;
*plot Current*Voltage=2/haxis=axis4 vaxis=axis3;
title "Time-ordered Residuals";
axis1 label=(a=90 'Residuals');
axis2 label=('Time Order');
*axis3 label=(a=90 'Current');
*axis4 label=('Voltage');
symbol1 v=dot h=.10 i=join ci=blue cv=blue;
*symbol2 v=dot h=.80 i=none cv=blue;
run; 
title;

proc sgpanel data=outa;
panelby trt / rows=1 columns=2;
series x=order y=resid/group=run;
rowaxis label="Residuals"; *This is applied to the rows statement, not the x axis;
colaxis label="Run Order"; *This is applied to the columns statement, not the y axis;
run;

/* Compound Symmetry */
proc mixed data=a;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
repeated / Subject=Run(Trt) type=cs;
run;

/* AR(1) model */
proc mixed data=a;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
repeated / Subject=Run(Trt) type=ar(1);
run;

/* GLIMMIX AR(1) */ 
proc glimmix data=a;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
*GLIMMIX uses _residual_ for R-side effects;
random _residual_/Subject=Run(Trt) type=ar(1);
run;

/* Toeplitz model; a banded Toeplitz model
could be added later (type=toep(q))*/
proc mixed data=a;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
repeated / Subject=Run(Trt) type=toep;
lsmeans trt voltage / pdiff adjust=tukey cl;
run;

/* GLIMMIX Toeplitz */ 
proc glimmix data=a plots=all;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
random _residual_/Subject=Run(Trt) type=toep;
lsmeans trt voltage / pdiff adjust=tukey cl;
run;

/* Unstructured model--this fails to converge */
proc mixed data=a;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
repeated / Subject=Run(Trt) type=un;
run;

/* GLIMMIX unstructured model--this fails to converge */
proc glimmix data=a plots=all;
class Trt Run Voltage;
model Current=Trt Voltage Trt*Voltage;
random _residual_/Subject=Run(Trt) type=un;
lsmeans trt voltage / pdiff adjust=tukey cl;
run;